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Abstract. To characterize the dynamical stability of an asteroid orbit, the calculation of 
its Lyapunov-Time Tl is a widely used mean in celestial mechanics. In the present work 
we investigated the effects of the used computer hardware and integration method on the 
outcome of such stability computations. We showed that for some asteroids the change of 
the employed numerical method can change the obtained Tl significantly. 

As a result of our investigations we introduced the computability index n as a measure 
of repeatability of such computations. 



1. Introduction 

Chaos indicators, like the Lyapunov exponent A, are widely used in celestial mechanics 
to characterize the dynamical behavior of bodies. The principal concept is to calculate 
the local exponential divergence of arbitrarily close initial conditions. Out of it the 
stability of their orbits can be determined. One might assume that a numerical 
calculation of A from the variational equations is straight forward. However, in the 
literature a lot of discrepancies between different studies dedicated to the same object 
can be found. 

In this work the dependency of such calculations on hardware and integration 
methods as a possible source of these differences is investigated. Being part of some 
bigger project we used here the MEGNO indicator to obtain A. A paper showing the 
sensitivity of Tl when being calculated from A directly is in preparation. 

2. Numerical Investigations 

2. 1. THE MEGNO CHAOS INDICATOR 

The MEGNO chaos indicator Y was defined for a dynamical system y = f(y) of 
dimension d by Cincotta et al. (2003) as 
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Here t denotes the overall integration time. 5 is the norm of the vector 8 in tangent 
space. Its time derivative 5 satisfies the variational equation 

The Lyapunov exponent A can be recovered from Y(t) by a least square fit, having 
the advantage of using the whole dynamical information contained in the integration 
interval. The Lyapunov time Tl is defined as the inverse of A. 

2. 2. SET-UP OF THE INVESTIGATION 

For this work a sample of real asteroids was used. These asteroids were chosen in a 
way, so that their Lyapunov times could be compared to already published values. 
Additionally, their Lyapunov exponent should be large enough to converge within the 
integration time span of 10 6 years to a finite value. So we decided for a sample of 3 
main belt asteroids: 

• Asteroid Helga (a = 3.633 AU and e = 0.074), whose Lyapunov time is given 
around 7000 years in all publications and 

• Asteroids 1981EY39 (a = 2.238 AU and e = 0.095) and 1994EL (a = 2.324 AU 
and e = 0.160), whose Lyapunov times were considered in Knezevic, Z. & 
Ninkovic, S. (2005) as difficult to determine. 

These asteroids were integrated together with a different number of planets of our 
solar system. All initial orbital elements were taken from the JPL HORIZONS system 
for the 01/01/2000. 

Table 1: Lyapunov times in 10 3 years for different integrations. 5 planets refer to 
the planets Mars to Neptune, while for 7 planets additionally Venus and Earth were 
integrated together with the asteroids. 

5 planets 7 planets 







Intel 


SGI 


Intel 


SGI 




ODEX 


6.69 


5.77 


5.87 


5.98 


Helga 


SABA 6 


10.51 


6.78 


6.04 


7.98 




SABA 8 


9.09 


7.50 


6.51 


7.45 




ODEX 


27.73 


29.59 


54.63 


116.80 


1981EY39 


SABA 6 


42.97 


52.62 


38.28 


42.56 




SABA 8 


40.57 


33.25 


73.07 


45.27 




ODEX 


57.95 


198.20 


49.31 


36.55 


1994EL 


SABA 6 


55.98 


74.14 


24.88 


53.37 




SABA 8 


43.59 


25.87 


34.66 


23.95 



The numerical integration of the equations of motion and their variational equa- 
tions was implemented in FORTRAN. As integrators the general purpose extra- 
polation scheme ODEX from Hairer et al. (1987) and the so called SABA symplectic 
scheme as described in Laskar & Robutel (2001) of order 6 and 8 were used. 
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For the calculation we had two different hardware platforms at our disposal: a 
SGI Origin3800 IRIX cluster with MIPS R120000 processor and a personal computer 
with an Intel Pentium 4 processor running under Linux. 

3. Results 

The results of all the different integrations are shown in Table [TJ One sees, that with 
the same program and the same initial conditions just e. g. by changing the used 
computer hardware, the result can change by a factor of 2 and more. 

The reasons for these differences arc inevitable numerical errors, like round-off and 
approximation errors, in the course of integration. These errors give rise to a slightly 
different trajectory, when changing the hardware architecture. And due to the long 
integration times to reach convergence of A, the same initial conditions could end up 
in very different regions of phase space at the end of the integration. 

By exploring the immediate neighborhood of the initial conditions, we found that 
the state of the nearby phase space determines the computability of Tl . If the phase 
space is homogeneous, slight changes in the initial conditions e. g. due to different 
rounding on different computers will still result in a similar Lyapunov time. For a 
coarse phase space (caused by overlapping resonances) these differences in rounding 
could result in very different trajectories of the asteroid and therefore different Tl as 
well. As an example the phase spaces in the plane of initial x and y coordinate for 
asteroid Helga and 1981EY39 are given in Figure [U 

Furthermore it was found that the general structure of the phase space is robust 
to a change of the used hardware platform as well as of the integration method. 
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Figure 1 : Color-coded stability map of the Lyapunov time Tl in the vicinity of asteroid 
1981EY39 (left) and Helga (right). The asteroids were integrated together with the 
planets Mars to Neptune for 10 6 years using OBEX. The initial cartesian x and 
y coordinates of these asteroids (point of origin) were changed in the last digit to 
produce a grid of 11 x 11 initial conditions of test particles, which are shown by black 
dots. To achieve a smooth representation the area between these points was linear 
interpolated. 
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Table 2: The computability index k for the asteroids Hclga, 1981EY39 and 1994EL. 





5 planets 


7 planets 




Intel 


SGI 


Intel 


SGI 


Hclga 


0.95 


0.97 


0.95 


0.97 


1981EY39 


0.95 


0.96 


0.72 


0.79 


1994EL 


0.87 


0.71 


0.93 


0.94 



4. The Computability Index k 

The information about the nearby phase space can be used to determine the repro- 
ducibility of a calculated Lyapunov time when changing the way of computation. To 
do so we generated a 5 x 5 grid of test particles by changing the last digit of the 
initial x and y coordinate of the asteroid. For these 25 points the corresponding Lya- 
punov times were calculated. Since the obtained results are not normally distributed, 
we used a method of robust statistics (the trimmed mean) to calculate the mean 
Lyapunov time Tr, and its standard deviation ct^t . 

As a measure of reproducibility we defined a computability index k as 

K = 1-JX (3) 

T L 

The closer this value is to 1, the more reliable is the result. The computability 
indexes for the investigated asteroids are shown in Table O 

5. Conclusions 

We showed that the used computer as well as the chosen integration method can 
influence the obtained Lyapunov time significantly. With the proposed computability 
index k the reliability of the stability time of a specific asteroid can be easily estimated. 

Finally we would like to remark that the commonly used tests to ensure the quality 
of an integration -the check of conservation of first integrals- give no indication of the 
just described problems. The error in energy stayed below 10~ n for all runs. 
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